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ABSTRACT 

We model the escape of ionizing radiation from high-redshift galaxies using high-resolution Adaptive 
Mesh Refinement A^-body+hydrodynamics simulations. Our simulations include time-dependent and spatially- 
resolved transfer of ionizing radiation in three dimensions, including effects of dust absorption. For galaxies 
of total mass M > 10" M0 and star formation rates SFR « 1-5 Moyr"', we find angular averaged escape 
fractions of 1 -3% over the entire redshift interval studied (3 < z < 9). In addition, we find that the escape 
fraction varies by more than an order of magnitude along different lines-of-sight within individual galaxies, 
from the largest values near galactic poles to the smallest along the galactic disk. The escape fraction declines 
steeply at lower masses and SFR. We show that the low values of escape fractions are due to a small fraction 
of young stars located just outside the edge of HI disk. This fraction, and hence the escape fraction, is progres- 
sively smaller in disks of smaller galaxies because their HI disks are thicker and more extended relative to the 
distribution of young stars compared to massive galaxies. Our results suggest that high-redshift galaxies are 
inefficient in releasing ionizing radiation produced by young stars into the intergalactic medium. We compare 
our predicted escape fraction of ionizing photons with previous results, and find a general agreement with both 
other simulation results and available direct detection measurements at z ~ 3. We also compare our simula- 
tions with a novel method to estimate the escape fraction in galaxies from the observed distribution of neutral 
hydrogen column densities along the lines of sights to long duration 7-ray bursts. Using this method we find 
escape fractions of the GRB host galaxies of 2 - 3%, consistent with our theoretical predictions. 

Subject headings: cosmology: theory - galaxies: dwarf - galaxies: evolution - galaxies: formation - stars: 
formation - methods: numerical 
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1. INTRODUCTION 

Star forming galaxies have a number of important effects 
on the surrounding intergalactic medium (IGM) and subse- 
quent gas accretion. The ionizing radiation from galaxies 
is thought to be responsible for the re-ionization of the uni- 
verse (e.g., Madau et al. 1999; Bolton et al. 2005), alterin 
the thermal state of the IGM (iGnedin & Hui 1998; GnediS |= 



20001: iRicotti et al."200tf,'Scha Ye et alFoOO^, McDonald et all 
200 It IFechner & Reimers 200' !), reducing gas accretion onto 
small dwarf galaxies (Efstat hioul 119921 : iThoul & Weinberg 
[199 6; Dijkstra et al. 2 004), and evaporating the existing gas 
in small hal os (Baik ana & Loeblll999l;ls"haviv & Deke]||2003l; 
IShapiro et al.i.2004) . 

The amount of radiation emitted by galaxies into the IGM 
depends not only on the abundance of hot, young stars, but 
also on the spatial distribution of absorbing gas and dust in 
individual galaxies and their immediate surroundings. The es- 
cape of ionizing radiation from galaxies is therefore the focus 
of a number of observational and theoretical studies. Never- 
theless, the "escape fraction", /esc, that characterizes the frac- 
tion of total ionizing radiation released into the IGM from 
individual galaxies remains poorly constrained. 

Recent empirical measurements of the escape fraction 
from normal galaxies in the local universe ([Leitherer et all 
[HHlHeckman et al."200lVBergvall et al.'2006) and at high 
redshifts (Giallongo et al. 2002; Fernandez-Soto et al. 2003; 
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IShaplev et all l2006h ^ have generally produced modest val- 
ues in the range of a few percent. In contrast, theoreti- 
cal studies of the escape of ionizing radiation from high- 
redshift galaxies have largely been inconclusive. Many of 
the previous studies have applied simplified analytic models 
(Dove & Shull 1994; Haiman & Loeb 1997; Dove et al .l200d: 
^cotti & Shull 2000; Wood & Loeb i2000e iClarke & Oeyl 
2002; Fujita et al. 2003) that predicted a wide range of val- 
ues for /esc. 

A more reliable estimate can be derived from self- 
consistent cosmological numerical simulations of galaxy for- 
mation that include three-dimensional radiative transfer and 
model the three-dimensional distribution of absorbing gas in 
and around individual galaxies. This approach is computa- 
tionally challenging and has been attempted o nly rec ently at a 
range of redshifts, from z ~ 20 ([Alv arez et al ll2006h to more 
modest redshifts z > 2 (Razoumov & Sommer- L arsenll2006i 
2007), where a comparison with the observations is possible. 

In this paper, we continue this line of work using fully 
self-consistent cosmological simulations that include radia- 
tive transfer and resolve the interstellar medium in modest- 
sized galaxies at high redshifts. Our work is similar but 
complementary to the work of Razoumov & Sommer-Larsei^ 
(l2006l l2007i) . In both approaches, self-consistent cosmo- 
logical sijnulations_ofjiojm 

while iRazoumov & Sommer-LarsenI (l2006i l2007h use the 
Smooth Particle Hydrodynamics (SPH) for modeling galax- 
ies with modest numbers of stellar particles (up to sev- 
eral thousand), we apply Adaptive Mesh Refinement (AMR) 
method in this study to reach a larger spatial dynamic range. 



^ Note that both 'Giallongo et al.I <2002l) and 'Sha plev et all <2006) quote 
larger values, up to 15%. However, the quoted numbers are for relative es- 
cape fractions between the 1500 and the Lyman limit. The absolute escape 
fractions are actually about 3%, as we discuss in more detail in ^ 
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iRazoumov & Sommer-LarseDl ( l2006ll200 7t) use the exact ray- 
tracing method for modehng the radiative transfer of ioniz- 
ing radiation from simulated galaxies without taking into ac- 
count effects of radiative transfer during simulation, while we 
use an approximate Opti cally Thin Variable Eddington Ten- 
sor (OTVET) method of Gnedin & Abel ( 2001) with radiative 
transfer calculations performed during the course of simula- 
tion self-consistently. 

Most importantly, in this work we explicitly include ab- 
sorption by dust, which allows a direct comparison with ob- 
servational estim ates of the escape fraction. G iven the dif- 
ferences between fRazoumov & Sommer-LarsenI (12006. 2007) 
and our approaches, a comparison of the results is useful for 
estimating the range of systematic numerical uncertainties in 
the problem, where a rigorous convergence study is not yet 
feasible. 

The paper is organized as follows. In §|2]we describe the 
simulations and numerical techniques used in our study. In 
§[3] we provide operational definition of the escape fraction 
and in § |4] we describe our model for dust absorption of the 
ionizing radiation. We present the results of the study in § |5] 
and discuss their implications in § |6] 

2. SIMULATION 

The simulation we use in this paper is described in detail by 
iTassi s et al. (2006, where it is labeled as FNEQ-RT). Here we 
briefly summarize its setup and numerical parameters. 

The simulation is performed using the Eulerian, gas dy- 
namics + A^-body Adaptive Refinement Tree (ART) code 
(iKravtsov et al.l 1 19971 l2002h . A large dynamic range is 
achieved through the use of Adaptive Mesh Refinement 
(AMR) in both the gas dynamics and gravity calculations. 
The calculation is started from a random realization of a Gaus- 
sian density field at z = 50 in a periodic box of 6/1"' Mpc in 
a flat ACDM model (Q,,, = 1-^^ = 0.3, = 0.043, h = 0.7, 
Hj = 1, and (7s = 0.9). A Lagrangian region corresponding to 
five virial radii of a Milky Way sized galaxy at z = is iden- 
tified and evolved with 2.6 x lO*" dark matter particles with 
masses of 9.2 X IO^/i-'Mq. 

The code employs a uniform 64^ grid to cover the entire 
computational box. The Lagrangian region is, however, al- 
ways unconditionally refined to the third refinement level, 
corresponding to the effective grid size of 512^'. As the matter 
distribution evolves, the code adaptively and recursively re- 
fines the mesh in high-density regions beyond the third level 
up to the maximum allowed 9th refinement level, which corre- 
sponds to the comoving spatial resolution of 260 pc, or phys- 
ical resolution of 65(50) pc at z = 3(4). 

The ART code computes metallicity-dependent, non- 
equilibrium gas cooling "on the fly" based on the abundances 
of five atomic species and molecular hydrogen. Star forma- 
tion and feedback (both radiative and thermal via supernova 
explosions and ste llar winds) are included, as described in 
ITassis et al.l (l2006l) . 

The self-consistent 3D radiative transfer of UV radiation 
from indivi dual stellar particles is followed with the OTVET 
algorithm (iGnedin & Abell l200ll) . The details of our im- 
plementation of the OTVET algorithm on adaptively refined 
meshes will be described elsewhere (Gnedin & Kravtsov 
2007, in preparation). A comparison of our radiative trans- 
fer scheme with several other existing numerical implemen- 
tations for the time-dependent a nd spatiall y-inhomogeneous 
radiative transfer is presented in llliev et al.l ([2006,) and in the 
foUow-up paper (Iliev et al. 2007, in preparation). 



In addition to the main progenitor of the Milky Way type 
galaxy, the fully refined Lagrangian region contains several 
dozens of smaller galaxies spanning a wide range of masses. 
We will use all these systems to estimate the escape fraction 
of ionizing radiation from star forming regions for galaxies of 
different masses and star formation rate (SFR). 

3. METHOD 

Because our simulation includes an approximate treatment 
of the radiative transfer self-consistently, the non-equilibrium 
abundances of all ions for each cell in the AMR grid hierarchy 
are available throughout the calculation. We can thus measure 
the opacity at a given frequency along a given direction from 
any point in the simulation volume to any other point. 

Since the opacity can only increase along a given line of 
sight, the escape fraction is, in general, a function of the dis- 
tance from the source. We choose to measure this distance in 
units of the virial radius of a given galaxy and we measure 
the escape fraction at 0.5, 1, and 2 virial radii from the center 
of the galaxy (defined as the peak of the dark matter density). 
Unless specified otherwise, hereafter we use the escape frac- 
tion at one virial radius as our fiducial working definition of 
the escape fraction. 

In computing the escape fraction, it is important to separate 
satellite galaxies from isolated ones. When a small satellite 
galaxy is located within a larger galaxy, the concept of virial 
radius and escape fraction becomes ambiguous. Therefore, 
we only consider isolated galaxies and only include sources 
that actually reside in (i.e. gravitationally bound to ) the main 
galaxy, thus excluding the satellites. In computing the contin- 
uum absorption by ionic species we include the detailed ve- 
locity structure of the galactic and infalling gas, although we 
find that ignoring the gas velocity field makes no measurable 
difference on the escape fraction. 

Operationally, the escape fraction /esc of a specific galaxy 
at a given redshift is then a function of two variables: the fre- 
quency and the direction of propagation of escaping radiation: 



/esc = fesc{v,9). 



(1) 



For studies of reionization of the universe and the Lyman- 
alpha forest, the most important quantity is the escape fraction 
of ionizing radiation. 



JLM = 7-, — 77 > 



(2) 



where al is the photo-ionization cross-section for the specie j 
(j = HI, He I, He II) and Si, is the spectrum of the sources of ra- 
diation. In our simulation, we include only stellar sources and 
compute the ionizing radiation spectra using the Starburst99 
package (Leitherer et al. 1999). We assume continuous star 
formation, 0.04 solar metallicity (typical of galaxies in our 
simulation), and a Salpeter initial mass function over a mass 
range fr om 1 to 100 Mp). Th e spectral shape is shown in Fig- 
ure 4 of iRicotti et alJ (l2002l) . Note that the Starburst99 spec- 
tral shape is computed for the unobscured stellar population, 
which is the appropriate spectral shape to use in equation 
which explicitly includes the /esc factor 

However, observationally the escape fraction of ionizing ra- 
diation is difficult to determine, since it requires measuring 
the whole ionizing spectrum. Instead, the escape fraction at 
the hydrogen ionization edge (Lyman limit), /esc(i^o), is usu- 
ally measured in observational studies. We present the rela- 
tionship between these two quantities in the appendix. 
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TABLE 1 

Parameters for the dust extinction law fits 



Terni 



A, (/im) 



Pi 



SMC dust model 



1 0.042 


185 


90 


2 


2 


2 0.08 


27 


15.5" 


4 


4 


3 0.22 


0.005 


-1.95 


2 


2 


4 9.7 


0.010 


-1.95 


2 


2 


5 18 


0.012 


-1.8 


2 


2 


6 25 


0.030 





2 


2 


7 0.067 


10 


1.9 


4 


15 


LMC dust model 


1 0.046 


90 


90 


2 


2 


2 0.08 


19 


21 


4.5 


4.5 


3 0.22 


0.023 


-1.95 


2 


2 


4 9.7 


0.005 


-1.95 


2 


2 


5 18 


0.006 


-1.8 


2 


2 


6 25 


0.02 





2 


2 


7 0.067 


10 


1.9 


4 


15 



' Values changed from IPe) 119921) are shown in bold, 
4. ABSORPTION BY DUST 

Incorporating dust absorption into the simulations for cal- 
culating the escape fraction is critical, because dust may con- 
tribu te significantly to the absor ption of ionizing radiation 
(e.g. lWeingartner & Drainell200ll) . Unfortunately, properties 
of dust in high-redshift galaxies are not known well. In partic- 
ular, dust absorption cross section depends on the dust com- 
position and grain size distribution, which is measured only 
in nearby galaxies. In our analysis, we adopt the dust ex- 
tinction curves for Large and Small Magellanic Clou ds (LMC 
and SMC respectively) froml Weingartner & Draind (1^01) as 
representative of high redshift galaxies, because low metal- 
licities of these galaxies are closer to typical metallicities of 
high-redshift galaxies in our simulation. 

A convenient parameterization for th e dust extinction law 
in LMC and SMC is given byiPe^ (Il992h . based on the eai'lier 
data. 



i=l 



where the fitting function / has the form 

f(x,a,b,p,q)= ■ 



(3) 



(4) 



We correct the fits o f Peil (Il992h using the newer data from 
IWeingartner & Draind ( 2001b bv adding the seventh term in 
equation (|3]l in order to account for the narrow and asymmet- 
ric shape of the FUV peak in the dust extinction law. We also 
change values for some parameters from those adopted by .Peii 
([T992|). The values for the parameters w e use are li sted in 
Table[Tl and we show in bold changes froml Peil(ll992h . 

The overall normalization for the dust cross section is de- 
termined by parameters (To^lmc = 5.6 x 10"^^ cm^ and (Tq.smc = 
1.1 X 10"^^ cm^. With these functional forms, the dust cross 
sections for both LMC and SMC fit the plotted curves of 
IWeingartner & Draind (1200 lb to a few percent accuracy. 

In order to limit the range of possible dust effects, we con- 
sider three extreme dust models: 1) a model with no dust at 
all, 2) a model in which we assume that the dust column den- 
sity scales with neutral (atomic and molecular) gas column 




100 
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Fig. 1. — The frequency dependence of the angular averaged escape 
fraction for the central galaxy at z = 3 for several dust models: no dust 
(solid line), SMC dust with instant sublimation (eq. (5), dotted line), LMC 
dust with instant sublimation (eq. (5), short-dashed Une), and SMC dust 
with no sublimation (eq. (6), long-dashed line). In addition, we also show 
with a short-long-dashed Hne the SMC no subhmation model with the dust 
cross section increased by a factor of 3, as a certain overestimate of the 
dust abso rption. The black poi nt with error-bars shows the extinction at 
1600 from lAdelbereer & Steidell <200a derived as Aisoo = 4.43 + 1.99/3 with 
/3 = -1.46 ±0.48). 

density. 



A^du 



Z 
Zo 



x(A?Hi + 2A?H,) 



(5) 



(where Z is the gas metallicity), and 3) a model where the dust 
column density scales with the total gas column density (both 
neutral and ionized). 



A^du.st = — xA^H, 
Zo 



(6) 



where Zq = 0.32 (0.2) is the g as-phase metallicity of the LMC 
(or SM C. IWeltv et al.lll997r[r9 99). Note that emission Hne 



studies (iPeimbert et alj l2000t TKeller & Wood 2006) usually 
indicate somewhat larger values for both metallicities than 
the absorption line metallicities that we adopt. If we adopted 
larger values for the LMC and SMC metallicities, the dust ef- 
fect on our measured escape fractions would be even smaller. 

Physically, equation (|5]l implies a complete instant sublima- 
tion of dust in the ionized gas, while equation ^ implies no 
dust sublimation at all. Of course, the truth lies somewhere 
in between, but these two extreme cases bracket the range of 
possibilities. 

Figure[T|shows the angular averaged escape fraction for the 
Milky Way progenitor galaxy at z = 3 as a function of fre- 
quency. We show several dust models, together with typ- 
ical reddening correction f or Lyman Break Galaxies from 
lAdelberger & Steidell (|2000|) . We also show the SMC no sub- 
limation dust model with dust cross section arbitrarily in- 
creased by a factor of 3. As we have mentioned above, the no 
sublimation model is likely to overestimate the effect of dust 
absorption. Since it is highly unlikely that uncertainties in our 
adopted value for the SMC metallicity and the dust extinction 
curve are as large as a factor of 3, this latter model serves as 
an absolute (and, likely, implausibly high) upper limit for the 
dust absorption. 

A rather unexpected feature of Figure [T] is that the escape 
fraction above the Lyman limit is almost independent of the 
dust absorption. In order to understand this phenomenon, we 
show in Figure |2] the distributions of escape fractions from 
all positions within the central galaxy in a fixed, randomly 
selected direction on the sky. The distributions have simi- 
lar shapes: they are weak functions of /esc for < /e,sc < 1 
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Fig. 2. — The luminosity weighted distribution of the escape fraction at 
the Lyman limit in a random direction on the sky from all positions within 
the central galaxy for the SMC instant sublimation dust model (our fiducial 
one). Three separate distributions are shown: absorption by gas only (dotted 
line), absorption by dust only (dashed Une), and total absorption (dust + gas, 
solid line). Notice that the gas-only distribution is essentially identical to the 
total distribution, because Zct'' ct*^ ' at the Lyman limit; so the dotted line 
coincides with the soHd one. The histograms are normalized to add up to \, 
so the spike at /esc ~ contains most of the data points. 

and exhibit a primary peak at /esc ~ and a secondary peak 
at /esc ~ 1 ■ Thus, the escaping radiation is produced not by 
sources in a semi-opaque medium, with each source attenu- 
ated by a similar amount, but by a small fraction of essen- 
tially unobscured sources. We discuss this point in detail in 
§|5]below (see Fig.[8]l. 

For a distribution like those in Figure |2l the effects of dust 
absorption can be easily understood, if we ignore the "translu- 
cent points" at < /esc < 1 and use a toy-model distribution 
for the gas-only /e.sc,Hi, 

P(/esc,Hl) = Ct5{\ -/e.sc,Hl) + (l -a)(5(/esc,Hl), (7) 

where a ^ 1 is constant and 5{x) is a Dirac delta function. 
The average escape fraction in this model is simply 

(/escHl) = / /esc,HlP(/esc.Hl)l^/esc,HI = Ct. (8) 

Jo 

If the dust absorption is included, then the escape fraction at 
each location is changed by a factor that depends on the dust 
opacity, 

/esc ~ /esc.HI^ — ^ 

At the locations where /cschi ~ 1 (thi ^1), the effect of dust 
is negligible (since ^ thi). At locations where Td > 1, the 
hydrogen opacity is already so large that no radiation escapes 
from this position, irrespective of how much dust is mixed 
with the gas. As the result, the escape fraction does not change 
at all. 

In reality, "translucent points" with < /esc < 1 are affected 
by dust, but their integrated contribution to the average escape 
fraction remains small. In the no sublimation dust model, the 
situation may be more complex, since there is dust absorption 
in the ionized gas. However, Figure [T| demonstrates that this 
is not a large effect. 

At frequencies below the Lyman limit the distribution of 
/esc from dust is similar to Figure |2l but there is no gas ab- 
sorption. Because dust absorption is weaker than gas absorp- 
tion, the average escape fraction at these frequencies is much 
larger and increases to unity as the frequency decreases down 
to infra-red. 

The important result of this section is that, while dust ab- 
sorption is the dominant effect for UV radiation below the 
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Fig. 3. — The escape fraction of HI ionizing radiation f^l as a function 
of direction, as visible from the center of the Milky Way type galaxy, at the 
virial radius of the galaxy at z = 5, 4, and 3. The sky is oriented in coordinates 
aligned with galactic gas disk such that the disk is horizontal in the middle of 
the plot. Note a different color scale for z = 5 image. 

Lyman Limit, it does not substantially affect the escape frac- 
tion of ionizing radiation in our model galaxies, as Figure [T] 
demonstrates. 

In the rest of this paper, we adopt the SMC instant subli- 
mation model as our fiducial dust model, but we note that our 
final results are not sensitive to this particular choice. 

5. RESULTS 

Figure [3] shows the escape fraction of hydrogen ionizing 
radiation f^J. (eq. |2]i as seen from the center of the central 
galaxy (the Milky Way progenitor) at three different redshifts. 
The celestial coordinates in Figure[3]are aligned with the prin- 
cipal axes of the galaxy ("galactic" coordinates), with the 
plane of the disk of the galaxy crossing the middle of the plot. 
Typically, the escape fraction is close to zero along the plane 
of the disk and approaches maximum values near the poles, 
although there are significant small-scale variations at all red- 
shifts that underscore the complex, perturbed nature of high 
redshift disk galaxies. At z = 5 the galaxy is experiencing a 
substantial major merger (and a lesser one at z = 3), so the an- 
gular distribution of escape fractions is more irregular at these 
epochs. 

Another interesting feature of the /esc distribution shown 
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Fig. 4. — The probability per unit log and per unit solid angle for having 
a specific value of the escape fraction of HI ionizing radiation for the central 
galaxy atz = 4 for 3 distances from the galaxy center (as measured in units of 
the galaxy virial radius). The shaded region around the middle curve shows 
the Poisson errors due to the finite number of pixels in the angular maps. 

in Figure [3] is a few small very opaque (dark blue) clouds of 
gas that block ionizing radiation at larger distances. These 
clouds can be counterparts of the Lyman Limit absorbers ob- 
served in the spectra of distant quasars. Previous studies, 
based on lower resolution simulations, have indicated that 
Lyman Limit systems te nd to cluster around large galaxies 
jKohler & Gned inll2007|). Our simulations confirm that such 
clouds do exist around high-redshift galaxies. 

We have checked that, as we integrate further in distance, 
more of Lyman Limit systems fall inside the radius of inte- 
gration, and the "sky", as seen from the center of the galaxy, 
appears progressively more opaque. The sky should become 
completely opaque at a distance of a few mean free paths 
for ionizing radiation. At z = 4 the mean fre e path for ioniz- 
ing radiation is about 85 co-moving Mpc ( iMiralda-Escudel 
120031) . much larger than the size of our computational box. 
Thus, we cannot actually reach the full opacity with our cur- 
rent simulation. However, the important point to make is that 
the mean free path is much larger than a virial radius of any 
galaxy at these redshifts, so the concept of the "escape frac- 
tion" is well defined and robust at intermediate redshifts. 

Figure|4]illustrates this point in a quantitative way. It shows 
the probability density for having a specific value of the es- 
cape fraction on a celestial sphere at different distances from 
the center of the main progenitor of the Milky Way-sized 
galaxy. As the distance form the galaxy increases, the tail 
of the distribution at low escape fractions grows, as Lyman 
limit systems cover a progressively larger fraction of the sky. 
However, at distance comparable to the virial radius of the 
galaxy the effect of Lyman limit systems is small, as can be 
seen from the plot: the fraction of the sky below f^^ = 0.003 
increases from 3% to only 6% as the distance increases from 
0.5/; viR to 2Rvi^. 

Poisson errors on the probability density shown in Figure]?] 
remain small for escape fractions well below lO""^, indicating 
that angular maps from Figure ]3] resolve the structure in the 
gas down to that level. Of course, the small-scale structure in 
the gas is limited by the finite resolution of our simulation. 

Figure ]5]presents the main result of this paper: the angular 
averaged escape fraction for hydrogen ionizing radiation as a 
function of galaxy mass or star formation rate at a range of 
redshifts. We only show results for the galaxies which are re- 
solved down to the maximum ninth level of refinement. This 
resolution criterion corresponds approximately to a minimum 
mass of 1O'"M0 (or > lO'* dark matter particles). A gen- 




SFR (Mg/yr) 

Fig. 5. — The angular averaged escape fraction of HI ionizing radiation 
as a function of galaxy mass (top) and star formation rate (bottom) at three 
redshifts (squares). The error-bars show the 10%-90% range for all possible 
directions. 




Fig. 6. — The evolution of the angular averaged escape fraction of HI 
ionizing radiation for the most galaxy (Mtot = 4 X 1O"M0 at z = 3). The 
shaded band shows the 10%-90% range for all possible directions. 

eral trend of increasing escape fraction with increasing galaxy 
mass and SFR is clearly observed: the escape fraction changes 
by approximately two orders of magnitude from « lO"'^ to 
« 0.02-0.05 for total masses between 10'" to 4 x 10" 
(or SFR from 0. 1 to « 7 Mq yr"' ). At the same time, the fig- 
ure also shows that there is little change with redshift, from 
z = 5 to z = 3, either in the values of the escape fraction or in 
the dominant trend with mass or SFR. 

Figure ]6] further illustrates the lack of redshift evolution in 
the escape fraction found in our simulations. The figure shows 
the average escape fraction of the most massive Milky Way 
progenitor in the simulation at different redshifts, along with 
the variation in escape fraction in different directions. While 
the escape fraction fluctuates with redshift, a spread in the 
escape fraction in different directions at a given redshift is 
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Fig. 7. — The mass dependence of the escape fraction of ionized radiation 
for tliree ionic species at z = 3 (similar to the top panel of Figure |3). We 
show only well-resolved galaxies and omit error-bars for clarity. The escape 
fractions are clumped at the level of 10"* to allow all points to be visible on 
the graph. 




Fig. 8. — The pseudo-color composite edge-on view of the main galaxy. 
The top half shows the gas density, while the bottom half depicts the neutral 
hydrogen fraction. Each image is composed of 3 snapshots from the sim- 
ulation taken lOMyr apart at z ~ 4 and combined as red, green, and blue 
channels of a color image. The predominance of purple (red and blue) color 
on top and green color on bottom indicates that the galactic disk oscillates 
globally on a time scale of about 20Myr. The oscillation pattern is much 
more chaotic in HI image outside of about 1 kpc. (This figure should be 
viewed in color) 

much larger than the variation in the average escape fraction 
with redshifts at z < 7. A similarly weak trend is found by 
most well resolved galaxies from our simulation. 

There is clearly a significant drop in the escape fraction at 
Hell ionization threshold. This is further illustrated in Figure 
|2l which shows all three escape fractions for ionizing radia- 
tion as a function of the total mass for model galaxies. We 
note that HI and He I are invariably close to each other, while 
Hell escape fractions are systematically much lower. This 
behavior is again completely expected, as only active galactic 
nuclei are thought to fully (i.e. doubly) ionize helium. 

We also study the mechanism for the escape of ionizing 
radiation from galactic disks in our simulations and for the 
predicted small values of /esc- As Figure |2] demonstrates, the 
ionizing radiation that escapes is emitted preferentially by un- 
obscured sources, with only a small fraction of all sources be- 
ing unobscured at any given time (rather than being emitted 



by the majority of sources, which are only partially obscured). 
In order to visualize the regions of a galaxy where these un- 
obscured sources reside. Figure [8] shows a pseudo-color com- 
posite image of three different snapshots from the simulation 
closely spaced in time. 

The edge-on disk of the main galaxy is shown in the to- 
tal gas density (top) and in the HI fraction (bottom). Three 
snapshots lOMyr apart at z 4 are overlaid together so that 
the first snapshot is shown in red, the second in green, and the 
last one in blue. If the gas distribution did not change between 
the snapshots, then the image would appear as an equal com- 
bination of red, green, and blue, i.e. as a pure grayscale. In 
reality, however, the upper surface of the disk appears purple 
(red plus blue), while the bottom surface is pure green. This 
means that the galactic disk oscillates with a period of about 
20Myr: the disk bended downwards between the first (red) 
snapshot and the second (green) snapshot, but returned to the 
first configuration in the third (blue) snapshot, so that the first 
and third snapshots look almost identical, blending into a sin- 
gle purple color. Oscillating modes with shorter wavelengths 
are also visible in the image. 

This behavior is clearly visible in the total gas density im- 
age and in the H I fraction image within the inner 0.5-0.7 kpc; 
the HI disk oscillates much more violently at larger radii. This 
is not surprising, because the oscillations of the HI edge are 
subject to the positive feedback. As the HI edge oscillates and 
uncovers ionizing sources, the ionizing intensity in a given 
point is likely to increase (similarly to how the ionizing inten- 
sity increases when isolated HII regions merge during reion- 
ization). 

Figure [8] thus shows why ionizing radiation of some young 
stars can leave the galaxy with virtually no attenuation. The 
cold HI disk, where young stars form is very thin (~ 100- 
200 pc). Some of the stars may have sufficient velocity to 
move outside the edge of the HI disk while they are still 
young, bright emitters of UV photons. For example, a star 
traveling at 10km s"' will move by w lOOpc in 10^ years. 
In addition, the outer edge of the HI disk is changing signif- 
icantly on the same time scale. These changes can expose 
some of the young stars and clear the way for the ionizing ra- 
diation to leave the system. The escaped ionizing radiation 
is thus mostly due to a small fraction of stars in a thin shell 
surrounding the HI disk. 

The relative constancy of the escape fraction with redshift 
or SFR that we find is then a simple consequence of the galac- 
tic disks oscillating by a comparable amount on lOMyr time 
scale. Note that this time scale is close to the local dynam- 
ical time scale of the dense high-redshift disks and lifetime 
of massive, UV-bright stars. Such oscillations are thus likely 
dynamical in origin: the disk can be perturbed by mergers 
and interactions with satellites, which are only weakly corre- 
lated with SFR or redshift (unless a disk goes into a starburst 
phase). 

This picture also explains why the escape fraction decreases 
steeply for lower mass galaxies. Gas distribution in dwarf 
galaxies is relatively more extended in comparison to the stel- 
lar distribution in our simulations, due to lower efficiency 
of star forrnation in dwarfs relative to massive galaxies (see 
iTassis et al]|2006l for detailed discussion). Stars are thus con- 
centrated towards central regions of the disk and a smaller 
fraction of them can reach the edge of HI disk or get exposed 
by fluctuations of its boundary. Note, however, that resolution 
of dwarf galaxies in our simulations is considerably worse 
than for the massive galaxies. This explanation and the trend 



7 



0.1 



10- 



10- 



10- 



10- 



I I I 

This work 

Razoumov &c Sommer — Larsen 2006 
Fernandez-Soto et al 2003 
Giallongo et al 2002 
Shapley et al 2006 




0.1 1 10 

SFR (M„/yr) 



10= 



Fig. 9. — A comparison of the angular averaged escape fraction at hydrogen 
ionization edge at z = 3 from various detemiinations. Filled red/gray upward 
triangle and circles show the obs ervational measurements from Shapley et al. 
{2006) and Femandez-Soto et al. (2003) respectively, presented as average 
values for their data samples. Two filled red/gray downward triangles show 
the upper limi ts from Giallongo et al. 1 2002) . Both Shapley et al. (2006) are 
IGiallongo et a l. ( 2002) points are corrected from the measured relative escape 
fraction to the absolute escape fraction, as explained in the text (open symbols 
show the original uncorrected measurements). The filled blue/black diamond 
is the simulation results of Razoumov & Sommer-Larsen (2006) for the sin- 
gle galaxy they report the star formation rate for. Filled blue/black squares 
with eiTor-bars show our results, similar to Figure[5] 

with mass will thus need to be verified in the future by higher 
resolution simulations. 

6. DISCUSSION AND CONCLUSIONS 

To summarize the current understanding of the escape frac- 
tion of photons at the ionizing threshold of hydrogen, we show 
in Figure|9]both simulation predictions and observational con- 
strai nts available at z = 3. The s imulation predictions are taken 
from lRazoumov & Sommer-L arsen (200 q) and o ur work. Ob- 
servational constraints are from Giallongo et al. (2002, upper 
limits only), Fernandez-Soto et al. (2003), and Shapley et al, 
(^006). The measurements of i Giallongo et al. (2002) are fo r 
individual galaxies wi th SFR taken from f ettini et ah (Il998h . 
The measurement of IShaplev et alj (I20()6l) is an average of 
14 star-forming galaxies at z ^ 3. The SFR is estimated 
based on the observed UV flux and corrected for dust ex- 
tinc tion for (E(B-y)) = O-H- Both IGiallongo et alj ( i20()l 
and IShaplev et al.l ( 120061) measurements are corrected from 
the relative measurements of the ratio /esc(912)//esc(1500) 
they report to the absolute measurement of /esc(912) by adopt- 
ing a value of /esc(1500) = 0.2 from Adelberger & SteideJ 
(|2000j). This correction factor is also cons istent with other 
observed estimates of reddening at 150 (Pettini et al.lll998l 
ISteidel et al.1 120011; ISh aplev et^ l2006l) . The measurements 
of iFernandez-Soto et al.i (i200ir are averaged over 14(13) 
high(low) luminosity galaxies found at 1.9 < z < 3.5 in the 
Hubble Deep Field (the second one is actually an upper 
limit). The star formation rates are conversio ns of the ob- 
served UV flux using the scaling relation from iMadau et al.l 
(IT998h . and are corrected for dust extinction assuming E(B- 
V) = 0.1 and the extinction law from Calze tti et aTl fcoOOl) . 
iFernandez-Soto et al.l (l2003h measure the absolute escape 
fraction, so no correction is applied to their points. 

There still exists a gap between simulations and observa- 




0.01 



1023 



Nh, (cm-=) 



Fig. 10. — A comparison of the cumulative distribution of HI column den- 
sities toward ionizing sources from our simulation at z = 3 (dashed line) with 
that measured in the host galaxies of long duration 7-ray bursts at z > 2. The 
black sohd line show s the observational data compiled by Vrees wiik et all 
120041) and IJakobsson et al. (2006), as well as the data for two additional 
sources from Shin et al. (2006) and Ruiz-Velasco et al. (2007). In two cases 
(GRB 060124 and GRB 060607) only upper limits have been measured. The 
dotted line is a power-law fit to the observational data for the range of HI 
column densities between lO"-^ cm^ and 10^'^ cm-. 

tions in the luminosity and SFR parameter space. While 
our simulations model modest-mass galaxies, observational 
measurements are reported mostly for the brightest, > L» 
galaxies. Nevertheless, the simulations are tantalizingly close 
to reaching a similar range of luminosity and SFR covered 
by the observations. Both the observations and our simula- 
tions indicate little (if any) trend in the absolute escape frac- 
tion with the galaxy mass or star formation rate, except for 
a rapid drop in /esc for smaller mass galaxies. Our simu- 
lations indicate that this characteristic scale corresponds to 
about SFR'^ 1 Mq yr"' or M tpt ~ 10" M,^, although an u p- 
per limit of /esc < 0.4% from lFernandez-Soto et al.l (l2003l) at 
SFR^ lOM0yr"' may indicate that simulations overestimate 
escape fractions in galaxies with SFR< lOMQyr"'. The level 
of the disagreement (if any) cannot yet b e accurately deduced 
from o ur simulations, and, formally, the iFernandez-Soto et al.l 
( I2OO3I) upper limit is fully consistent with our results. Note 
also that Si ana et al.l (l2007l) in their recent study of the Lyman 
continuum escape fractions in galaxies at z ~ 1 also find small 
values (/esc ^ 0.1), which is consistent with weak evolution 
of escape fractions we find in our simulations. The general 
agreement between our results and observations is definitely 
encouraging and may indicate that relative distribution of the 
young stars and UV absorbing gas is modeled faithfully in our 
high-resolution AMR simulations. 

An interesting alternative way of constraining the escape 
fraction observationally, independent of the observed UV 
light, is offered by the observed distribution of HI column 
densities in the host galaxies of long duration 7-ray bursts 
(GRBs; Chen, Prochaska, & Gnedin 2007, in preparation). 
Assuming that GRBs sample the distribution of young stars in 
an unbiased fashion, the fraction of young stars located in ion- 
ized regions with small or negligible attenuation of UV radi- 
ation should correspond to the fraction of GRBs with HI col- 
umn densities lower than some threshold value corresponding 
to the transition between neutral and ionized hydrogen. 

Figure [Tol shows the observed distribution of HI column 
densities in spectra of GRBs with confirmed redshifts of z > 2, 
compared to the corresponding distribution in our simulation 
at z = 3. The observed distribution is remarkably close to 
a power-law for 10''''^cm^ < A^hi < 10^'-'' cm^. Extrapolat- 
ing this power-law to lower column densities of « (5 - 8) x 
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lO'^ cm^ (or Til = 3-5, appropriate for the boundary between 
the ionized and neutral gas^) impHes that about 2- 3% of all 
GRBs are located in ionized regions in which the UV radia- 
tion of their progenitors should have escaped freely into the 
IGM, so that 

Asc « /grb(A^hi < (5-8) X lO'^cm^) = 0.02-0.03. 

This value is remarkably close to the escape fractions of 
« 1-3% measured in our simulations. This test also pro- 
vides strong support for our interpretation of why the escape 
fractions are so low. 

Note that the H I column density distribution toward ioniz- 
ing sources (i.e. column densities from source locations to 
one virial radius, as explained in §3) from our fiducial simu- 
lation at z = 3 is reasonably close to the observed distribution 
of A^Hi from GRBs for A^hi < 3 x 10^' cm^. At higher column 
density we do not expect a good agreement between our sim- 
ulation and the data, because the simulation does not incorpo- 
rate the physics of formation and self-shielding of molecular 
hydrogen correctly, which becomes important in that regime. 
Thus, the simulation is expected to overpredict A^hi column 
densities, because in reality some of this HI is in molecular 
form. 

While higher resolution simulations and deeper observa- 
tions are needed to bridge the remaining small gap in the 
probed star formation regime between the data and the theory, 
our results suggest that average escape fractions from bright 
galaxies at intermediate redshifts do not depend strongly on 
galaxy properties or redshift. This is interesting because this 
implies that porosity of the interstellar medium in galaxies, 
and hence its opacity to ionizing radiation, does not dramat- 
ically change at higher star formation rates, as expected in 

* The specific location of tlie ionization edge of a galaxy depends on a 
number of factors and does not correspond to a particular value of th. A 
value of TLL = 1, which corresponds to the suppression of ionizing back- 
ground by only a factor of 3, is certainly not enough to make the gas sub- 
stantially neutral. A value of th = 10 is, on the other hand, enough to make 



some of the theoretical models (e.g., IClarke & Oe"vll2002h . 
Note that observed escape fractions at higher star formation 
rates are consistent with being simply an extrapolation of the 
trend seen in simulations at lower rates. If this is indeed the 
case, this implies absence of a well-defined critical SFR above 
which the escape fraction sharply increases to unity. 

At the same time, the low values of escape fractions found 
in our simulations suggest that high-redshift galaxies are quite 
inefficient in emitting ionizing radiation produced by their 
young stars into the intergalactic medium, with escape frac- 
tions decreasing sharply for galaxies of Mtot ^ 10" Mq (or 
SFR< \Mq yr"'). This conclusion potentially has important 
implications for the contribution from normal galaxies to the 
early reionization of the universe and for the relative role of 
galaxies and quasars in keeping the universe ionized at inter- 
mediate redshifts. 



We thank Douglas Rudd for constructive and useful com- 
ments on the draft of this paper This work was supported in 
part by the DOE and the NASA grant NAG 5-10842 at Fer- 
milab, by the NSF grants AST-0206216, AST-0239759, and 
AST-0507596, and by the Kavli Institute for Cosmological 
Physics at the University of Chicago. Supercomputer simula- 
tions were run on the IBM P690 array at the National Center 
for Supercomputing Applications and San Diego Supercom- 
puting Center (under grant AST-020018N) and the Sanssouci 
computing cluster at the Astrophysikalisches Institut Pots- 
dam. This work made extensive use of the NASA Astro- 
physics Dat a System, arXiv . org preprint server, and the 
HEALPix^ (lG6rskietalJl2005h package. 

the gas even at the mean cosmic density largely neutral. Thus, a value of tll 
around 3 to 5 should correspond to the approximate location for the galaxy 
ionization edge. 

'^ |http : / /healpix . jpl . nasa ■ gov| 



APPENDIX 

ESCAPE FRACTION FOR IONIZING RADIATION AND THE ESCAPE FRACTION AT THE IONIZATION EDGE 

In order to facilitate comparison between the theoretical and observational results, in Figure[TT]we show a relationship between 
the theoretically relevant escape fraction of ionizing radiation (i.e. an integral quantity, as defined in eq.|2l since it determines the 
ionization state of the IGM) and the escape fraction at the hydrogen ionization edge (Lyman limit), which is usually measured in 
observational studies. We find a tight correlation between the two quantities in the form 

/»i=1.25/e.c(j^0), (Al) 

which is helpful for relating observationally measured and theoretically relevant quantities. We use this relation in Fig. |9] when 
comparing observational and theoretical values on the same plot. 
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